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ABSTRACT 

We employ hydro dynamic simulations to study the effects of high shock compression 
ratios, as expected for fast shocks with efhcient particle acceleration, on the convective 
instability of driven waves in supernova remnants. We find that the instability itself 
does not depend significantly on the compression ratio, a, with the growth rates and 
the width of the mixing region at saturation being comparable for the range of ratios 
we studied; 4 < o" < 21. However, because the width of the interaction region between 
the forward and reverse shocks can shrink significantly with increasing a, we find that 
convective instabilities can reach all the way to the forward shock front if compression 
ratios are high enough. Thus, if supernova blast waves accelerate particles efficiently, 
we expect the forward shock to be perturbed with small amplitude, small wavelength 
bumps, and to find clumps and filaments of dense ejecta material in the vicinity of the 
shock. In addition and in contrast to situations where o" < 4, any enhancement of the 
radial magnetic field from Rayleigh- Taylor instabilities will also extend all the way to 
the shock front and this may help explain the slight dominance of radial fields long seen 
in polarization measurements of young remnants like Tycho. 

Subject headings: hydrodynamics - instabilities - shock waves - ISM: supernova remnants: 
cosmic rays 



1. INTRODUCTION 

The presence of convective or Rayleigh- 
Taylor (R-T) instabilities in young, super- 
nova remnants (SNRs) has been argued for 
a long time (e.g., Guh 1973; Shirkey 1978; 
Dickel et al. 1989) and examined in detail 
with numerical simulations of idealized mod- 
els (e.g., Chevalier et al. 1992; Jun & Nor- 
man 1996a; Wang & Chevalier 2000). These 
instabilities are expected to arise in young. 



ejecta-dominated SNRs where a shell of dense, 
shocked ejecta is gradually decelerated by 
lower-density shocked circumstellar material. 
This situation is subject to the familiar R-T 
instability, which results in fingers of dense 
ejecta gas protruding into, and mixing with 
the shocked circumstellar material. These in- 
stabilities are relevant to remnants from both 
Type la and Type II explosions, and are ex- 
pected to exist as long as there is a reverse 
shock present in the remnant. 



These instabilities may be important for 
mixing ejecta out to, and possibly past, the 
forward shock, and for modifying the mor- 
phology of an otherwise spherical outer blast 
wave and/or reverse shock. One of the best 
studied young SNRs, Gas A exhibits ejecta 
material out to and beyond the nominal ra- 
dius of the forward shock, as evidenced by 
optical knots of enriched gas (Fcscn & Gun- 
dcrson 1996), as well as spectral analysis of X- 
ray knots and filaments (Hughes et al. 2000a). 
The outer blastwave of Gas A, as seen in X-ray 
emission, is roughly spherical, but possesses 
many bumps and protrusions that may be 
the result of the convective instability. Other 
examples of irregular blastwaves include the 
more extreme cases of SN1986J and SN1993J. 
VLBI observations of these SNRs (Bartel et 
al. 1991; Bietenholz, Bartel, & Rupen 2001) 
show shells of radio emission with large ra- 
dial protrusions, suggesting that very young 
SNRs can exhibit large deviations from spher- 
ical symmetry. 

Another characteristic of young SNRs that 
may originate with the convective instability 
associated with the contact discontinuity is 
the slight dominance of a radial magnetic field 
as deduced from radio polarization measure- 
ments (Dickel et al. 1991). Jun &; Norman 
(1996b) suggested that the polarization ob- 
servations could be explained by convective 
instabilities dragging out an ambient mag- 
netic field as the R-T fingers pushed from the 
contact discontinuity out toward the forward 
shock. 

Despite the belief that convective insta- 
bilities alone should play an important role 
in mixing ejecta out to the forward shock 
and generating the observed radial magnetic 
fields, attempts to demonstrate this have 
failed. Ghevalier et al. (1992) used a linear 
perturbation analysis and two-dimensional 
hydrodynamic simulations to study this in- 
stability, and found that the mixing region 



did not reach the forward shock. This result 
was confirmed by Jun & Norman (1996a) and 
found to hold true even in the presence of 
rapid cooling of the shocked ejecta (Gheva- 
lier & Blondin 1995). Jun & Norman (1996b) 
showed that this limited mixing region was in- 
consistent with observations, for it produced 
a SNR with a two-shell structure when viewed 
in radio-synchrotron emission: an inner shell 
of strong emission and a dominant radial mag- 
netic field, and an outer shell of weaker emis- 
sion with a dominant tangential field. Jun k. 
Jones (1999) extended this work by including 
a simplified scheme for injecting and acceler- 
ating test-particle electrons and calculated the 
radio synchrotron emission self-consistently in 
the turbulent fields. Their work emphasized 
the role of the reverse shock in accelerating 
electrons and producing radio emission. 

One possible solution to these problems 
is the introduction of small-scale structure 
into the density of either the circumstellar 
medium (GSM) or the supernova ejecta. Jun, 
Jones, &: Norman (1996) found that including 
small cloudlets in the surrounding GSM could 
enhance the growth of the Rayleigh- Taylor 
fingers enough so they reached all the way 
to the shock front. Alternatively, Blondin, 
Borkowski, k, Reynolds (2001) examined the 
effects of including low-density bubbles in the 
supernova ejecta. They also found an en- 
hanced turbulence in the region between the 
forward and reverse shocks, with some clumps 
of ejecta reaching out to the forward shock. 

All of this previous work concerning R-T 
instabilities and driven waves has assumed the 
supernova ejecta and GSM can be treated as 
an ideal gas, typically modeled with an adi- 
abatic index, 7 = 5/3.^ However, these as- 
sumptions may not be appropriate for young 
SNRs, where the strong shock waves are ex- 

^For discussions of other types of instabilities in shocks 
accelerating cosmic rays see, for example, Drury & 
Falle (1986) and McClements et al. (1996). 
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pected to be efficient accelerators of energetic 
particles (i.e., cosmic rays). ^ It has long been 
believed that SNRs are the primary sources 
of galactic cosmic rays below ~ lO-*^^ eV (e.g., 
Axford 1981). However, for this to be the 
case, the acceleration mechanism must be ef- 
ficient and place on the order of 10% or more 
of the total ejecta kinetic energy into relativis- 
tic particles (e.g., Blandford &: Eichler 1987; 
Drury, Markiewicz, & Volk 1989; Dorfi 2000). 
Radio observations have long been proof that 
SNRs produce GeV electrons and recent evi- 
dence of X-ray synchrotron emission in young 
SNRs suggests they can produce TeV elec- 
trons as well (e.g., Koyama et al. 1995; Mas- 
tichiadis & de Jager 1996; Allen et al. 1997; 
Tanimori et al. 1998; Reynolds &: Keohane 
1999). While direct proof of ion accelera- 
tion in SNRs remains elusive, indirect evi- 
dence exists (e.g., Elhson, Slane Sz Gaensler 
2001) and there is a growing body of evidence 
linking the efficient production of cosmic-ray 
ions by the forward and reverse shocks with 
the thermal properties of the shock-heated X- 
ray emitting gas (e.g., Dorfi &; Bohringer 1993; 
Decourchelle, Ellison, & Ballet 2000; Hughes, 

^In the heliosphere, where shocks are observed directly 
with spacecraft, there is clear evidence that the quasi- 
parallel earth bow shock (with a sonic Mach number, 
Mgo < 10) can place 10-30% of the solar wind kinetic 
energy flux into superthermal particles (e.g., Ellison et 
al. 1990). Quasi-perpendicular interplanetary shocks 
(IPSs) also accelerate ambient particles but generally 
with lower efficiencies due to their lower Mach num- 
bers (generally MgQ ^ 3 for IPSs) (e.g., Baring et 
al. 1997). In a few CEises, however, IPSs with higher 
than average Mach numbers have been observed to 
accelerate particles with efficiencies comparable to the 
quasi-parallel Earth bow shock (Eichler 1981; Tera- 
sawa et al. 1999). Hybrid plasma simulations showing 
direct injection and acceleration of ambient particles 
are reasonably consistent with these observations (e.g., 
Scholer et al. 1992), as are convection-diffusion mod- 
els (e.g., Kang & Jones 1997). Roughly, quasi-parallel 
and quasi-perpendicular shocks are those with an an- 
gle between the upstream magnetic field and the shock 
normal of less than 45° and greater than 45° , respec- 
tively. 



Rakowski, & Decourchelle 2000b). 

One of the important structural effects of 
efficient particle acceleration is a shock com- 
pression ratio greater than 4 (e.g., Ellison Sz 
Eichler 1984; Berezhko & Ellison 1999). For 
an ideal gas in the absence of particle acceler- 
ation, a high Mach number shock produces a 
compression ratio. 



For the nominal value of 7 = 5/3, a = 4. 
However, when efficient particle acceleration 
occurs, relativistic particles can be created 
and substantially contribute to the pressure, 
and superthermal particles can escape from 
the shock. These two effects combine to pro- 
duce a shocked plasma which acts to a large 
extent like a gas with an effective adiabatic 
index, 7eff < 5/3, allowing arbitrarily large 
compression ratios. In fact, it has been shown 
(Berezhko &; Ellison 1999) that, when injec- 
tion of particles from the background occurs 
easily, nonlinear diffusive shock acceleration 
yields 

r 1.3 if 1<M|o<Mao ; 

(7 ~ < 

[1.5 MX if 1<Mao<M|o . 

(2) 

Here, Mgo (Mao) is the upstream sonic 
(Alfven) Mach number (see also Kazanas & 
Ellison 1986; Malkov 1997). Furthermore, 
the increase in compression ratio is accompa- 
nied by a decrease in the post-shock temper- 
ature which links the efficient production of 
superthermal particles with the thermal prop- 
erties of the shock-heated, X-ray emitting gas. 

Of particular importance for the work re- 
ported here, the increased compression ratio 
also implies that the region between the for- 
ward and reverse shocks in young SNRs will 
be significantly narrower and denser when ac- 
celeration is efficient compared to the case 
where little energy is placed in relativistic par- 
ticles. The increased density and narrower 
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spatial extent of the interaction region leads 
to large density gradients suggesting that R-T 
instabilities may have faster growth rates and 
larger amplitudes at saturation than is com- 
monly assumed (Decourchelle, Ellison, & Bal- 
let 2000). The dependence of the compression 
ratio on the sonic and Alfven Mach numbers 
implies that these effects will be most pro- 
nounced in SNRs having high shock speeds 
and relatively low magnetic fields. 

In this paper we explore the possibility that 
the familiar convective instability in young 
SNRs may be altered by high compression ra- 
tios. We use a well-tested hydrodynamic sim- 
ulation of evolving SNRs, only replacing the 
standard 7 = 5/3 with an effective adiabatic 
index, j^g < 5/3. A softer equation of state (a 
low 7eff) causes the shocked plasma to be more 
compressible, yielding shock compression ra- 
tios C7 > 4, consistent with those expected 
when particle acceleration is efficient. 

Somewhat unexpectedly, we find that the 
convective instability growth rate and mixing 
length are almost independent of the compres- 
sion ratio. However, since the region between 
the forward and reverse shocks narrows as a 
increases, the likelihood that Rayleigh- Taylor 
fingers reach and perturb the forward shock 
is strongly correlated with a and, therefore, 
with the efficiency of cosmic-ray production. 
If a is large enough, the convective instabili- 
ties do reach and perturb the forward shock 
making it likely that clumps and filaments of 
dense ejecta material will be found there. 

In section 2 we discuss the impact of a low 
value of 7off on a one-dimensional dynamical 
model of young SNRs. In section 3 we present 
hydrodynamic simulations in one, two, and 
three dimensions, and in section 4 we discuss 
the approximations and implications of our 
numerical results. 



2. A DYNAMICAL MODEL 

The dynamical evolution of young SNRs 
is determined by the interaction of the stel- 
lar ejecta, one to several solar masses of ma- 
terial thrown out by the SN explosion, with 
the surrounding circumstellar medium. In the 
simplest case, this interaction is characterized 
by a two shock structure: a forward shock 
driven ahead of the ejecta into the CSM, and 
a reverse shock that decelerates the super- 
sonically expanding ejecta in the vicinity of 
the interaction region. To model this cjccta- 
dominated evolution, we begin with the self- 
similar driven wave (SSDW) model (Chevalier 
1982a), although we have found results simi- 
lar to those reported here for the exponential 
ejecta model (Dwarkadas & Chevalier 1998). 

Chevalier (1982a) provided spherically 
symmetric, self-similar solutions to the struc- 
ture of cjccta-dominated SNRs under the as- 
sumption that the mass density in both the 
ejecta and the CSM could be described by 
power laws. Since the steep power law den- 
sity profiles for the ejecta lead to an infinite 
ejecta mass, we follow Chevalier and insti- 
tute a cutoff in the ejecta density profile at 
small radii. Therefore, at radii less than some 
critical value, Tc = Vct, the ejecta density is 
taken to be a spatially constant plateau, while 
beyond this radius it follows a power law, i.e.. 



Pct~ 



for r > Vct 
for r < Vct 



(3) 



where t is time, the constant pc is given by 
5n-25 



27rn 



-EsnV' 



(4) 



the velocity at the intersection of the ejecta 
density plateau and the power law is given by 

1/2 

Vr = I — I , (5) 



3n-9 Me- 

and n is a constant index. These quantities, 
Pc and Vc, are related to the constant g in 
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Chevalier (1982a) by p^v^ = 5". Here, E.^ 
and Moj are the kinetic energy and mass of 
the ejecta, respectively, and it's clear that fix- 
ing the plateau allows us to relate the self- 
similar parameters to physical quantities nor- 
mally associated with SNRs. Furthermore, 
we restrict our examples to times such that 
Tc is less than the radius of the reverse shock 
to avoid complications from multiply reflected 
shocks etc., and to make direct comparisons 
between our numerical results and self-similar 
solutions straightforward. 

The density of the ambient circumstellar 

medium is also described by a power law in 
radius in the Chevalier (1982a,b) solution, i.e., 



Pa = qr 



(6) 



where s = corresponds to a uniform den- 
sity (in which case q is the ambient density) 
and s = 2 represents the density profile of a 
steady-state stellar wind from the supernova 
progenitor of speed, Vyj, and mass loss rate, 
dM/dt, such that q = {dM / dt) / {Air v^) . 

Thus, provided that the edge of the ejecta 
plateau region has not yet reached the reverse 
shock, the radius and expansion velocity of 
the forward shock are: 



Ri 



Vi 



PcVc \ {n-z)/{n-s) 



n-S\ R 



n — s J t 



, (7) 

(8) 



and the fiuid variables in the shocked region 
immediately behind the forward shock are: 



Pi 



Ul 



Pi 



7eflF+ 1 



qRi' , 



Teff - 1 

2 / n-3 \ Ri 
7eflF-|- 1 \n- sj t 



2q 



7eflF-|- 1 \n- s 



n-SV Rl 



(9) 
(10) 

, (11) 



where p is the density, u is the fluid speed, p 
is the pressure, and the subscript 1 indicates 



the forward shock. Note that all of the de- 
pendence on the adiabatic index in eq. (7) is 
contained in the constant b, which is related 
to parameters defined in Chevalier (1982a) 
through the expression b = {Ri/Rc)A^^^'^~^\ 
The calculated values of the constant A and 
the ratio of the forward shock radius to the ra- 
dius of the contact discontinuity for 7 = 5/3 
for various n and s are given in Table 1 of 
Chevalier (1982a). 

We note that in an attempt to predict 
the 7-ray flux from pion-decay. Chevalier 
(1983) considered two-fluid, self-similar solu- 
tions where one fluid was a thermal gas with 
7 = 5/3 and the other was a relativistic gas 
(i.e., cosmic rays) with 7 = 4/3. The effective 
adiabatic index was determined by 



7eff - 



5 + 3w 



(12) 



where w was the (constant) fraction of total 
pressure made up by relativistic particles at 
the shock. Chevalier did not address the pos- 
sible effects 7cff would have on R-T instabili- 
ties nor attempt to model the effects of par- 
ticle escape during acceleration and assumed 
4/3 < 7ofi' < 5/3 with a maximum cr = 7 in 
his calculations. Chevalier did note, however, 
that cosmic rays diffusing from the shocked 
gas might produce a dense shell with a > 7 
similar to that which occurs in the later ra- 
diative phase. 

The structure of a SSDW is illustrated in 
Figure 1, which shows the radial profiles of 
density, pressure, and velocity for values of 
n = 7 and s = 2. The interaction region is 
bounded by the forward and reverse shocks, 
which for this example have compression ra- 
tios (7 = 4 corresponding to a high Mach 
number adiabatic shock in an ideal fluid with 
7 = 5/3. The pressure is monotonically in- 
creasing from the reverse shock to the forward 
shock, as all of the shocked gas is gradually 
decelerating. 
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Fig. 1. — Internal structure of a SSDW with 

n = 7, s = 2, and 7 = 5/3 propagating to 
the right. The contact discontinuity separates 
the shocked ejecta on the left and the shocked 
CSM on the right. The radius is scaled to the 
position of the forward shock and the density, 
pressure, and velocity are scaled to their val- 
ues just inside the forward shock. 

An important characteristic of these solu- 
tions for s = 2 is the convectively unsta- 
ble region ahead of the contact discontinuity 
(Chevalier et al. 1992). This region is marked 
by opposing signs of the pressure and density 
gradients. The positive pressure gradient is 
gradually decelerating the shocked ejecta and 
shocked CSM, but the negative density gradi- 
ent implies that the denser gas is being decel- 
erated by lower density gas; a situation that is 
subject to the familiar R-T instability. A sim- 
ilar situation holds for s = solutions, but in 
this case the opposing gradients are behind 
the contact discontinuity rather than in front. 

While the idealized self-similar hydrody- 
namical model has been very successful in 
describing general aspects of SNR evolution, 
it is clearly incomplete in some important 
ways. For instance, the hydrodynamic model 
ignores any effects from a magnetic field, the 
one-dimensional approximation eliminates the 
ability to describe irregularities in the ambi- 



ent conditions or Rayleigh- Taylor generated 
instabilities, and particle acceleration at the 
collisionless forward and reverse shocks is 
neglected. Here we focus on how efficient 
cosmic-ray production influences R-T insta- 
bilities by setting < 5/3. We find that 
lowering 7eff causes the blast wave to ex- 
pand slightly slower, but the radius of the 
forward shock still follows the self-similar re- 
sult Ri oc t("~3)/(n--s)_ More importantly, 
lowering 7eff causes the width of the region 
between the forward and reverse shocks to 
shrink considerably and alters the internal 
structure of the self-similar driven wave. 

We list some of the parameters describing 
SSDWs for a range of 7eff in Table 1 and, for 
7 = 5/3, these values are identical to those in 
Table 1 of Chevalier (1982a). These param- 
eters, as well as the radial profiles shown in 
Figure 2, were generated following the proce- 
dure outlined in Chevalier (1982a). Note in 
particular that the value of b, which contains 
the 7eff dependence, changes by less than 30% 
as the compression ratio is varied from 4 to 
21. Thus the velocity of the blast wave is not 
strongly affected by the value of the compres- 
sion ratio. In contrast, between o" = 4 and 21, 
the width of the interaction region shrinks by 
a factor of 4 in the case of n = 7 and s = 0. 

The radial profiles of SSDWs for different 
values of are shown in Figure 2. While 
the increase in the compression ratio and de- 
crease in the width of the interaction region 
are clearly evident, there is also a change in 
the density gradient of the shocked CSM. As 
gas moves from the forward shock to the con- 
tact discontinuity it compresses/expands adi- 
abatically to match the relatively constant 
pressure of the interaction region. A lower 7cff 
leads to a slower drop in pressure due to radial 
expansion, so the gas does not need to be com- 
pressed as much as with a higher j^s- The re- 
sult is a more positive density gradient in the 
shocked CSM for a smaller 7eff. From Figure 
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Fig. 2. — Radial profile of the gas density in 
SSDWs with n = 7 and s = 2 for values of %fi 
as labeled. Note the change in sign of the den- 
sity gradient behind the forward shock. Here 
and in Figure 3, the density is scaled to its 
value immediately ahead of the forward shock. 



2 we see that this can even lead to a change 
in the sign of the density gradient, which may 
in turn affect the stability properties of the 
waves. Nonetheless, we expect the low 7off 
driven wave to remain unstable, as there is 
still a shell of dense shocked ejecta being de- 
celerated by lower-density shocked CSM. 

3. HYDRODYNAMIC SIMULATIONS 

To examine the effects of different com- 
pression ratios on the convective stability of 
driven waves, we use the Virginia Hydrody- 
namics (VH-1) time-dependent hydrodynamic 
numerical code in one, two, and three dimen- 
sions. All of the simulations were computed in 
spherical geometry on a numerical grid with 
500 zones in each dimension. This resolution 
is comparable to the highest resolution sim- 
ulations of Chevalier et al. (1992). Based on 
their results, we expect that a higher resolu- 
tion simulation would show more small scale 
mixing, but the large-scale dynamics of the 
problem would not change. The angular ex- 



tent of the simulation domain was chosen to 
produce roughly square numerical zones, i.e., 
RAO « AR. Periodic boundary conditions 
were applied in both the 9 and (f) directions 
for the multidimensional simulations. The ra- 
dial boundary conditions were set to match 
the relevant power laws. The numerical grid 
was expanded to follow the forward blast wave 
so the evolution could be tracked for many ex- 
pansion times. 

We ran one-dimensional (ID) and two- 
dimensional (2D) simulations with the four 
combinations of s = 0, 2 and n = 7, 12. For 
each (n, s) pair we ran four simulations with 
7(,ff = 5/3,4/3,1.2, and 1.1. In the following 
discussions we will focus on the set of simu- 
lations with n = 7 and s = 2 as a specific 
example, but the results are qualitatively the 
same for all of these parameter choices. 

3.1. 1-D SIMULATIONS 

We ran ID simulations for all the cases 
listed in Table 1 as a check on our abil- 
ity to numerically recreate the analytic, self- 
similar solutions. The ID hydrodynamic sim- 
ulations were initialized with appropriate den- 
sity power laws, and allowed to evolve until a 
self-similar state was reached. The gas pres- 
sure in the ejecta and CSM was kept suffi- 
ciently low to ensure that the Mach numbers 
of the forward and reverse shocks were always 
greater than 100. 

In Figure 3 we show the results of our 
most extreme case with n = 12 and s = 
and 7 = 1.1 overlayed on the semi-analytic 
solution obtained by numerically integrating 
an ordinary differential equation (Chevalier 
1982a). This particular SSDW is very thin 
and possesses steep gradients, making it a dif- 
ficult test case for the hydrodynamic code. 
Nonetheless, our simulation does a reasonably 
good job of matching the analytic solution; 
note in particular the relative sharpness of 
the shock fronts and contact discontinuity. In 
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Fig. 3. — The density profile of a driven wave 

for the relatively extreme case with n = 12, 
s = 0, and 7cfr = 1.1. The solid line is the 
self-similar analytic solution and the triangles 
are from the ID hydrodynamic simulation. 

the less extreme cases, the match between the 
simulation and the analytic solution is better 
than shown in Figure 3. However, the very 
sharp density peak at the contact discontinu- 
ity in the s = 2 solutions, and the density 
trough in the s = solutions, are typically 
smoothed out over a few zones by this hydro- 
dynamic method. 

3.2. 2-D SIMULATIONS 

The multi-dimensional simulations were 
initiated with the analytic SSDW solutions 
normalized to i?i = 1, Vi = 1, and a preshock 
density of unity. The simulation thus be- 
gins at an initial time of t = (n — 3)/(n — s) 
in these normalized units. The initial solu- 
tion is mapped onto a spherical grid with a 
radial span roughly twice the width of the 
interaction region. The convective instability 
was seeded by adding acoustic noise to the 
shocked interaction region. These simulations 
were evolved long enough for the convective 
instability to reach saturation, such that the 
mixing region is in a quasi-self-similar state 
with the growth of R-T fingers matched by 



their destruction due to shearing and advec- 
tion. This typically required ~ 6-8 doubling 
times of the blast wave radius, and beyond 
this time, the results are qualitatively inde- 
pendent of the seed noise. 

In all of the cases listed in Table 1, the 
layer of shocked ejecta was found to be un- 
stable with an evolution similar to that seen 
in previous studies (i.e., Chevalier et al. 1992; 
Jun &; Norman 1996a). As an example, we 
show the density for driven waves with n = 7, 
s = 2, and different values of 7eff in Figure 
4. The results for other combinations of n 
and s arc similar. These images are taken 
at the end of the simulations, once the insta- 
bility has reached saturation. In addition to 
shrinking the interaction region, a lower value 
of 7eff produces a narrower, denser shell of 
shocked ejecta and sharper, denser R-T fin- 
gers. Aside from the narrower interaction re- 
gion, these results for small j^s are similar to 
those found by Chevalier & Blondin (1995) 
for driven waves with radiative cooling in the 
shocked ejecta. In this latter case the ejecta 
behaved as an ideal gas with 7 1, while 
the CSM was modeled as an ideal gas with 
7 = 5/3. 
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Table 1: Numerically derived parameters describing spherically symmetric driven waves. The sub- 
script 1 (2) indicates values immediately behind the forward (reverse) shock. 
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Fig. 4. — Convective instabilities in driven waves with n = 7 and s = 2 and values of 7off as marked. 
The color scale represents gas density, scaled to the density of the CSM immediately ahead of the 
forward shock. The thin shell of dense, shocked ejecta is deformed into narrow fingers, characteristic 
of the R-T instability. 
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We can obtain a rough estimate of the in- 
stability growth rate by following the turbu- 
lent energy density in the intershock region. 
Rather than attempt to subtract off the bulk 
radial velocity to find the contribution to the 
turbulent velocity, and to avoid confusion in 
comparing 2D and 3D simulations, we con- 
sider only a single angular component of the 
turbulent motion. Furthermore, to remove 
the effects of radial expansion, deceleration, 
and shock compression, we normalize this tur- 
bulent energy to the bulk kinetic energy den- 
sity associated with the shock front. We thus 
define the ratio 



(13) 



where the integration is only over the volume 
of shocked gas. The evolution of this turbu- 
lent energy parameter is shown in Figure 5, 
from which one can see that while the satu- 
rated levels of turbulent energy density differ 
by about a factor of two, the growth rates 
are not noticeably affected by changes in 7eff, 
despite the fact that the radial profile from 
the spherical solution has changed dramati- 
cally (e.g., the density gradient has changed 
sign). 

The saturation of the instability, at least in 
terms of the width of the mixing region, is also 
relatively unaffected by changes in 7eff. In Fig- 
ure 6 we plot an angle-averaged radial profile 
of the ejecta mass fraction. In a spherically- 
symmetric model this would be a step func- 
tion from unity dropping to zero at the con- 
tact discontinuity. To provide an easy com- 
parison, we have normalized the radii such 
that the reverse shock is located at the same 
radius for all models. From this we can see 
that the length of the R-T fingers is compa- 
rable for all four values of 7eff. 

While the instability itself is insensitive to 
changes in the compression ratio, the overall 
width of the interaction region decreases with 
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Fig. 5. — Growth of the convective instabil- 
ity as measured by the normalized turbulent 
energy density, the interaction region. 

These curves are from 2D simulations with 
n = 7, s = 2, and 7eff = 5/3 (solid), 4/3 (long 
dash), 1.2 (short dash), and 1.1 (dotted). 
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Fig. 6. — The width of the mixing region, as 
illustrated with the mass fraction of ejecta, is 
relatively unaffected by changes in 7eff. For 
this comparison we have normalized radii in 
each model to the reverse shock rather than 
the forward shock, so we can directly compare 
the radial extent of the mixing region. The 
line styles are as in Figure 5. 
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increasing compression ratio. This has the 
important effect that, for small 7cff, the R-T 
fingers can reach as far as the forward shock 
front. To measure the maximum extent of the 
R-T fingers, we identify an average mass frac- 
tion of 0.001 as the leading edge of the mixing 
region. Since the radial gradient of the mass 
fraction is quite steep at these small values, 
the exact choice of a cutoff value does not af- 
fect the results. The time dependence of the 
maximum radial extent (relative to the for- 
ward shock) of the R-T fingers is illustrated 
in Figure 7. In the run with the highest com- 
pression ratio (7eff = 1-1), the R-T fingers 
quickly reach and pass the average shock ra- 
dius. We note that when the curves in Figure 
5 become approximately flat, the system is in 
a near self-similar state and the ratio of the 
forward shock radius to the contact disconti- 
nuity radius remains approximately constant. 
If we allowed the system to evolve to the point 
where the reverse shock entered the plateau 
region in the ejecta density (which nominally 
happens when the mass of swept up CSM is 
comparable to the ejecta mass), the ratio of 
the forward shock radius to the contact dis- 
continuity radius would start to increase and 
the R-T fingers would drop behind the for- 
ward shock. 

Once the fingers of dense, shocked ejecta 
reach the forward shock, they can begin to 
distort the outer blast wave by pushing small 
regions out ahead of the average shock radius. 
The resulting bumps in the outer shock can be 
seen in Figure 4 for the simulation with 7cff = 
1.1. The shock front is nearly spherical for 
large values of 7efi, while values close to unity 
produce significant deviations from spherical 
symmetry. 

We note that even in the most dramatic 
case of n = 12, s = 0, and 7eff = 1.1, the de- 
viation from a spherical outer shock remains 
relatively small; the penetrating fingers arc 
not able to dramatically alter the forward 
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Fig. 7. — The maximum radial extent (mea- 
sured relative to the average radius of the 
outer shock) of the R-T fingers of dense 
shocked ejecta is shown as a function of time 
for n = 7 and s = 2. For the smallest value of 
7cff, the fingers reach all the way to the shock 
front. For the other values of 7off the R-T fin- 
gers reach a maximum length after t ~ 1000. 



shock. This can be seen in Figure 7, where 
the maximum radial extent of the fingers ap- 
pears capped at ~ R\. Since the instability 
is relatively unaffected by the value of 7eff, 
wc expect Figure 7 to show a series of par- 
allel curves as the fingers grow until satura- 
tion occurs (at several 100 times the initial 
simulation time, although this depends on the 
magnitude of the initial perturbation). If the 
fingers reach the forward shock before satu- 
ration, as in the 7eff = 1.1 case, they cannot 
push out much beyond the shock and their 
growth stalls. This limited ability to affect 
the shock front is due to the strong shearing 
flow created when the shock front is distorted. 
If a clump of dense ejecta has sufficient radial 
momentum to push out the forward shock in 
a local protrusion, the deformed shock front 
generates a substantial tangential post shock 
flow around the ejecta clump. In all of our 
simulations, this shear flow quickly disrupts 
the clump through the Kelvin-Helmholtz lu- 
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stability, with the remnants of the clump 
quickly advected back into the interaction re- 
gion. Thus we never see a local protrusion 
of the shock front stick out more than a few 
percent of the blast wave radius. 

3.3. 3-D SIMULATIONS 

We repeated two simulations in 3D to check 
that these results are not dramatically af- 
fected by the assumption of axisymmetry, us- 
ing n = 12 and s = for jcif = 5/3 and 1.1. 
In comparing the 2D and 3D simulations, we 
stress that the numerical code remained vir- 
tually identical. The only changes involved 
using a 3D grid and repeating the hydrody- 
namic updates in a third direction. 

For the more familiar case of %{[ = 5/3 
our results are quite similar to previous 3D 
simulations of driven waves (e.g. , Jun & Nor- 
man 1996b), showing only minor differences 
between 2D and 3D simulations. This similar- 
ity is exhibited in Figure 8. The only signifi- 
cant difference visible in this Figure is the in- 
creased amount of small-scale structure in the 
3D simulation, although this difference does 
not appear to affect any of the global proper- 
ties of the driven wave. A quantitative com- 
parison of these two simulations finds that the 
growth and saturation of the turbulent energy 
is virtually identical, the radial extent of the 
fingers is comparable (the 3D fingers reached 
slightly further than the 2D fingers), the for- 
ward shock remains spherical and the width 
of the interaction region grows slightly in both 
cases as the instability approaches saturation. 

The situation is somewhat more compli- 
cated for the case of high shock compression, 
as shown in Figure 9. Here one sees much 
more small-scale structure in the 3D simula- 
tion, both within the interaction region and 
in the forward and reverse shocks. This more 
dramatic difference between 2D and 3D in the 
high-compression case is consistent with the 
fact that the R-T fingers do not reach the 



shock front when the shock compression is 
low, but they do reach - and perturb - the 
forward shock when the compression is large. 
For both values of 7eff we see more small-scale 
structure in the R-T fingers in 3D, but in the 
high-compression case this small-scale struc- 
ture can modify the forward shock to produce 
many small wavelength, but large amplitude, 
perturbations. 

A consequence of the more oblique shocks 
created by the deformation of the forward 
shock in the 3D simulation is an overall drop 
in the average compression ratio and a corre- 
sponding increase in the width of the inter- 
action region. Despite these differences, our 
statistical measures of the convective instabil- 
ity remain relatively unaffected by the dimen- 
sionality of the simulation. Furthermore, our 
primary conclusion from the 2D simulations 
still holds in 3D; the R-T fingers are able to 
deform the forward shock when the compres- 
sion ratio is high. 
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Fig. 8. — A comparison of the shock structure in a 2D axisymmetric simulation (left), with a planar 
slice from a 3D simulation (right). The shading depicts gas density, scaled to the density of the 
preshock gas. Black represents the high density of the shocked ejecta and white the low density of 
the unshocked CSM. These runs used n = 12, s = 0, and 7ofr = 5/3. 
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4. DISCUSSION 

If young supernova remnant shocks accel- 
erate cosmic rays with high efficiency, the ac- 
celeration process can cause these high Mach 
number shocks to have compression ratios 
considerably greater than four. We have in- 
vestigated the effects of this high compres- 
sion on the hydrodynamic stability of young 
remnants with a simple and direct method; 
we perform hydrodynamic simulations with 
the effective adiabatic index, 7eff, set to val- 
ues less than 5/3, yielding compression ra- 
tios, a ~ (7eff + l)/(7eff - 1) > 4. While 
this fluid approach is clearly an approxima- 
tion to the real situation where the collision- 
less shocks accelerate particles to relativistic 
energies, we believe it accounts for the major 
effects of efficient particle acceleration on the 
overall dynamics and, most importantly, on 
the Rayleigh- Taylor instabilities that develop 
in the interaction region between the forward 
and reverse shocks. 

Somewhat surprisingly, we find that a 
large compression ratio has little effect on the 
growth of R-T instabilities. However, since a 
large a dramatically shrinks the width of the 
interaction region between the forward and re- 
verse shocks, several important effects emerge 
from of our high-compression simulations: 

(1) In contrast to shocks with a ~ 4, R-T 
fingers reach closer to the forward shock front, 
and if a is large enough, ejecta material can be 
found at, and even slightly ahead of the aver- 
age shock radius. This can occur fairly quickly 
after the explosion (depending on how quickly 
a reverse shock forms and the magnitude of 
any inhomogcneitics seeding this instability) 
and the fingers should stay near the forward 
shock front throughout the time the reverse 
shock is in the power law portion of the ejecta 
density profile. In general this should apply 
to SNRs from an age of only months up to 
several 1000 years. 



(2) The forward shock is perturbed on 
short wavelengths and with relatively small 
amplitudes. In addition to slightly altering 
the spherical shape of the shock front, this 
will lead to some small spread in the post- 
shock temperatures. Note also that the over- 
all temperature of the shocked gas will be sub- 
stantially lower if particle acceleration is effi- 
cient and compression ratios are large than if 
little acceleration occurs (e.g., Ellison et al. 
2000). This effect provides a coupling be- 
tween the shock morphology and thermal X- 
ray emission (Decourchelle, Ellison, &; Bal- 
let 2000; Hughes, Rakowski, & Decourchelle 
2000b). 

(3) If the mixing region reaches the for- 
ward shock, the morphology of the magnetic 
field as seen through radio synchrotron polar- 
ization observations may be affected. Jun & 
Norman (1996b) followed the evolution of the 
ambient magnetic field in driven wave simula- 
tions, looking for an explanation for the ori- 
gin of observed polarization in young SNRs in 
the elongation of the field by R-T instabilities. 
Our work offers a ready explanation for why 
this polarization can extend all the way to the 
forward shock. 

While we do not include magnetic fields 
in the hydrodynamic simulations we perform 
here, we note that any ambient magnetic fields 
will also be compressed at the shock front. 
In addition, the increased turbulence associ- 
ated with the high compression driven waves 
may be expected to further amplify the mag- 
netic fields and show stronger radio emission 
and possibly more intense TeV emission from 
inverse-Compton and/or pion-decay than the 
more quiescent driven waves found for ^y^s = 
5/3. It can be expected that the effects we 
see here from increased compression will add 
to those reported by Jun &; Jones (1999). 

The SNR IE 0102.2-7219 in the Small 
Magellanic Cloud may represent a young rem- 
nant for which this model of a high com- 
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prcssion driven wave is applicable. Recently, 
Hughes, Rakowski, & Decourchelle (2000b) 
have determined postshock electron temper- 
atures in SNR IE 0102.2-7219 using Chan- 
dra X-Ray Observatory observations. Using 
the measured forward shock speed of Vps ~ 
6000 km s~^, they find that the electron tem- 
perature of 0.4-1 keV is at least 2.5 times 
lower than can be explained with standard 
(i.e., test-particle) shock heating even if only 
Coulomb electron heating occurred. They 
conclude that the forward shock is placing at 
least 50% of the shock kinetic energy flux in 
cosmic ray ions. 
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Fig. 10. — The radial profile of density {p/nip) 
for a model of IE 0102.2-7219 with effi- 
cient particle acceleration assuming "f^s =1-2 
(dashed curve), and one with test-particle ac- 
celeration assuming 7efi = 5/3 (solid curve), 
both at an age of tgnr — 700 yrs. 

If one applies a SSDW model to the ob- 
served radius and expansion velocity of IE 
0102.2-7219 (e.g., Gaetz et al. 2000; Hughes, 
Rakowski, & Decourchelle 2000b) assuming 
n = 9 and s = 0, one finds an age of 700 years, 
independent of the compression ratio. Choos- 
ing reasonable parameters of = 10^^ erg 
and Moj = IM©, one finds an ambient density 
of Pa ~ 0.1 mp cm~^for £7 = 4. If instead we 
assume a large compression ratio of cr = 11, 
the derived ambient density is lower by only 



a factor of 2.3. Despite the relative similar- 
ity in fitting the observed SNR parameters, 
these models, shown in Figure 10 are dramat- 
ically different in terms of the position of the 
reverse shock and the postshock temperature. 
If the shock compression is high, the reverse 
shock will be significantly closer to the for- 
ward shock, as shown in Figure 10. This is 
not consistent with X-ray observations which 
show the presence of shocked cjccta at a radius 
of only ~ 70% of the forward shock (Gaetz et 
al. 2000). However, if the ejecta mass is rela- 
tively small, this SNR may be evolving away 
from the SSDW phase, with the reverse shock 
now propagating in toward the center of the 
SNR. 

With these values and one additional pa- 
rameter, the ambient magnetic field strength 
Bq, we can estimate the particle accelera- 
tion using the model of Berezhko & Ellison 
(1999). We find that with 5o ~ 3 /xG, the 
forward shock obtains a ~ 11 (consistent with 
7eff = 1.2) with a shocked proton tempera- 
ture of Tp2 ~ 4x 10^ K, i.e., 10 times lower 
than obtained in a test-particle shock with the 
same parameters (for the acceleration calcu- 
lation we assume a constant Vps = 6000 km 
s~^ for ~ 700 yr). The low postshock electron 
temperature deduced from X-ray observations 
is now easily explained if the electron temper- 
ature is ~ 1/3 Tp2.'^ Furthermore, the forward 
shock accelerates electrons to ^ 70 TeV in 
~ 700 yr, consistent with radio synchrotron 
emission. At the current age, Mso — 3700 
and Mao — 150, and the acceleration is ex- 
tremely efficient, i.e., more than 80% of the 
kinetic ram energy flux is placed in relativis- 
tic ions. While this model is not unique, it 
does show that efficient particle acceleration 
is consistent with reasonable supernova and 
ambient medium parameters, as well as with 
deduced values of radius, speed, age, and elec- 

*If we had used a lower 7efj, we could have obtained a 
consistent fit with a smaller Tp2. 
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tron temperature. 

While there are clearly differences between 
the actual situation in SNRs where some frac- 
tion of swept-up material is shock accelerated 
to relativistic energies, and simply lowering 
7ofi', we believe our results are qualitatively 
correct. The important differences include the 
fact that lowering 7eff implies that the effects 
of relativistic particle pressure and particle 
loss occur everywhere rather than just in the 
shocked gas and in the precursor regions in 
front of the forward and reverse shocks. We 
note that if the finite size of the shock is lim- 
iting the maximum particle energy, there will 
be a significant fraction of the total pressure 
in high energy particles with upstream diffu- 
sion lengths on the order of 1/10 of the shock 
radius. This difference in pressure distribu- 
tion will have some effect on the shock evo- 
lution since, instead of a uniform 7eff, the ac- 
tual remnant has a gas with a soft equation 
of state pushing one with a harder equation of 
state. While we have not tested this difference 
in detail, we expect the effects we found will 
actually be enhanced if 7eflF varies spatially. 

Another difference comes about because we 
use eq. (1) with a constant 7cff to determine 
a instead of eqs. (2). With (1), there is little 
variation in a as the SNR evolves as long as 
Mso 3> 1, while eqs. (2) can give a much larger 
variation depending on the parameters. ^ In 
any case, since in all probability a S> 4 during 
the time we consider before the reverse shock 
enters the ejecta plateau region, any differ- 
ences resulting from a varying with time are 
likely to be small and go in the direction of 
increasing the effects we report. 

"'We note that it is not always the case that lowering 
the Mach numbers causes a to decrease. Equations (2) 
apply when injection into the shock acceleration mech- 
anism is efficient. If injection is weak enough, high 
Mach number, test-particle solutions with cr ~ 4 can 
result (see Berezhko & Ellison 1999, for a full discus- 
sion) . 



Finally, since the shrinking of the interac- 
tion region between the forward and reverse 
shocks depends totally on the pressure in that 
region, effects other than particle acceleration 
which influence the pressure may be impor- 
tant. The most likely effect which we have 
neglected comes from the compression of the 
magnetic field, B. Since the magnetic pres- 
sure scales as B^, increasing the compression 
ratio could produce magnetic pressures large 
enough to prevent the interaction region from 
becoming narrow enough to allow the R-T 
fingers to reach the forward shock. This, of 
course, will depend on the Alfven Mach num- 
ber and the angle the upstream field makes 
with the shock. The pressure effects of the 
magnetic field will be offset somewhat by the 
fact that, for a given compression ratio, the 
shocked thermal pressure is considerably less 
in a shock undergoing efficient acceleration 
compared to one with a low 7eff and no ac- 
celeration. 
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